%-------------------------------------------------------------------------------

% This file is part of code_saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2025 EDF S.A.
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA 02110-1301, USA.

%-------------------------------------------------------------------------------

\programme{cfener}

\hypertarget{cfener}{}

\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Fonction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Pour les notations et l'algorithme dans son ensemble,
on se reportera à \fort{cfbase}.

Après masse (acoustique) et quantité de mouvement,
on considère un dernier pas fractionnaire (de $t^{**}$ à $t^{***}$)
au cours duquel seule varie l'énergie totale $E = \rho e$.

\begin{equation}\label{Cfbl_Cfener_eq_energie_cfener}
\left\{\begin{array}{l}
\rho^{***}=\rho^{**}=\rho^{n+1}\\
\\
\vect{Q}^{***}=\vect{Q}^{**}=\vect{Q}^{n+1}\\
\\
\displaystyle\frac{\partial \rho e}{\partial t}
+ \divs\left( \vect{Q}_{ac} \left(e+\displaystyle\frac{P}{\rho}\right) \right)
= \rho\vect{f}_v\cdot\vect{u}
+ \divs(\tens{\Sigma}^v \vect{u})
- \divs{\,\vect{\Phi}_s} + \rho\Phi_v
\end{array}\right.
\end{equation}

Pour conserver la positivité de l'énergie, il est indispensable ici,
comme pour les scalaires, d'utiliser le flux de masse convectif acoustique
$\vect{Q}_{ac}^{n+1}$ compatible avec l'équation de la masse.
De plus, pour obtenir des propriétés de positivité sur les scalaires,
un schéma upwind pour le terme convectif doit être utilisé
(mais les termes sources introduisent des contraintes supplémentaires
qui peuvent être prépondérantes et gênantes).

\vspace{0.5cm}

à la fin de cette étape, on actualise éventuellement
(mais par défaut non)
une deuxième et dernière fois la pression
en utilisant la loi d'état pour obtenir la pression finale~:
\begin{equation}
\displaystyle P^{n+1}=P(\rho^{n+1},\varepsilon^{n+1})
\end{equation}

See the \doxygenfile{cfener_8f90.html}{programmers reference of the dedicated subroutine} for further details.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discrétisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%---------------------------------
\subsection*{Discrétisation en temps}
%---------------------------------

La modélisation des flux de chaleur choisie jusqu'à présent est de la
forme $-\divs(\,\vect{\Phi}_s) = \divs(\lambda \gradv{T})$.

Pour faire apparaître un terme diffusif stabilisant dans la
matrice de résolution, on cherche à exprimer le flux diffusif de chaleur
($-\divs(\,\vect{\Phi}_s)$)
en fonction de la variable résolue (l'énergie totale).

Avec $\varepsilon_{sup}(P,\rho)$
dépendant de la loi d'état, on exprime l'énergie totale de la façon suivante~:
\begin{equation}
e = \varepsilon + \frac{1}{2} u^2
= (C_v T + \varepsilon_{sup}) + \frac{1}{2} u^2
\end{equation}

En supposant $C_v$ constant\footnote{Pour $C_v$ non constant, les
développements restent à faire~: on pourra se
reporter à  P. Mathon, F. Archambeau, J.-M. Hérard : "Implantation d'un
algorithme compressible dans \CS", HI-83/03/016/A}, on a alors~:
\begin{equation}\label{Cfbl_Cfener_eq_flux_thermique_cfener}
-\divs(\,\vect{\Phi}_s)
= \divs(K \gradv(e - \frac{1}{2} u^2 - \varepsilon_{sup}))\qquad
\text{avec } K=\lambda / C_v
\end{equation}

Lorsqu'un modèle de turbulence est activé, on conserve la même forme
de modélisation pour les flux thermiques et $K$ intègre alors
la diffusivité turbulente. On pourra se reporter à
la documentation de \fort{cfxtcl} à ce sujet.

Avec la formulation~(\ref{Cfbl_Cfener_eq_flux_thermique_cfener}),
on peut donc impliciter le terme en $\gradv{e}$.

\bigskip
De plus, puisque la vitesse a
déjà été résolue, on implicite également le terme en
$\gradv{\frac{1}{2} u^2}$. L'exposant $n+\frac{1}{2}$ de $\varepsilon_{sup}$
indique que l'implicitation de ce terme est partielle (elle dépend de la forme
de la loi d'état).

Par ailleurs, on implicite le terme de convection, le terme de puissance
des forces volumiques, éventuellement le terme de puissance des forces de
pression (suivant la valeur de \var{IGRDPP}, on utilise la prédiction de
pression obtenue après résolution de l'équation portant sur la masse
volumique ou bien la pression du pas de temps précédent)
et le terme de puissance des forces visqueuses. On implicite le terme de puissance
volumique en utilisant $\rho^{n+1}$.

\bigskip
On obtient alors l'équation discrète portant sur $e$~:
\begin{equation}\label{Cfbl_Cfener_eq_energie_totale_cfener}
\begin{array}{l}
\displaystyle\frac{(\rho e)^{n+1} - (\rho e)^n}{\Delta t^n}
+ \divs(\vect{Q}_{ac}^{n+1} e^{n+1}) - \divs(K^n \gradv{e^{n+1}})
= \rho^{n+1} \vect{f}_v \cdot \vect{u}^{n+1}
- \divs(\vect{Q}_{ac}^{n+1} \displaystyle\frac{\widetilde{P}}{\rho^{n+1}} )\\
\text{\ \ \ \ }+ \divs((\tens{\Sigma}^v)^{n+1} \vect{u}^{n+1})
- \divs(K^n \gradv(\frac{1}{2} (u^2)^{n+1}
+ \varepsilon_{sup}^{n+\frac{1}{2}}))
+ \rho^{n+1}\Phi_v\\
\end{array}
\end{equation}
avec $\widetilde{P}=P^{Pred}\text{ ou }P^n$ suivant la valeur de \var{IGRDPP}
($P^n$ par défaut).

En pratique, dans \CS, on résout cette équation en faisant apparaître à
gauche l'écart $e^{n+1} - e^n$. Pour cela, on écrit la dérivée
en temps discrète sous la forme suivante~:

\begin{equation}
\begin{array}{ll}
\displaystyle
\frac{(\rho e)^{n+1} - (\rho e)^n}{\Delta t^n}
& =
\displaystyle
\frac{\rho^{n+1}\, e^{n+1} - \rho^n\, e^n}{\Delta t^n}\\
& =
\displaystyle
\frac{\rho^{n}\, e^{n+1} - \rho^n\, e^n}{\Delta t^n}+
\frac{\rho^{n+1}\, e^{n+1} - \rho^n\, e^{n+1}}{\Delta t^n}\\
& =
\displaystyle
\frac{\rho^{n}}{\Delta t^n}\left(e^{n+1} - e^n\right)+
e^{n+1}\frac{\rho^{n+1} - \rho^n}{\Delta t^n}
\end{array}
\end{equation}

et l'on utilise l'équation de la masse discrète pour écrire~:
\begin{equation}
\displaystyle
\frac{(\rho e)^{n+1} - (\rho e)^n}{\Delta t^n}
=
\frac{\rho^{n}}{\Delta t^n}\left(e^{n+1} - e^n\right)-
e^{n+1}\dive\,\vect{Q}_{ac}^{n+1}
\end{equation}



%---------------------------------
\subsection*{Discrétisation en espace}
%---------------------------------


%.................................
\subsubsection*{Introduction}
%.................................

On intègre l'équation (\ref{Cfbl_Cfener_eq_energie_totale_cfener})
sur la cellule $i$ de volume $\Omega_i$ et l'on procède comme
pour l'équation de la masse et de la quantité de mouvement.

On obtient alors l'équation discrète
suivante~:
\begin{equation}\label{Cfbl_Cfener_eq_energie_totale_discrete_cfener}
\begin{array}{l}
\displaystyle\frac{\Omega_i}{\Delta t^n}
(\rho_i^{n+1} e_i^{n+1}-\rho_i^n e_i^n)
+ \displaystyle\sum\limits_{j\in V(i)}
\left(e^{n+1} \vect{Q}_{ac}^{n+1}\right)_{ij} \cdot \vect{S}_{ij}
- \displaystyle\sum\limits_{j\in V(i)}
\left(K^n\gradv(e^{n+1})\right)_{ij}\cdot\vect{S}_{ij}\\
\\
\text{\ \ \ \ } = \Omega_i\rho_i^{n+1} {\vect{f}_v}_i \cdot \vect{u}_i^{n+1}
- \displaystyle\sum\limits_{j\in V(i)}
\left(\displaystyle\frac{P^{Pred}}{\rho^{n+1}}\
\vect{Q}_{ac}^{n+1}\right)_{ij} \cdot \vect{S}_{ij}
+ \displaystyle\sum\limits_{j\in V(i)}
\left((\tens{\Sigma}^v)^{n+1} \vect{u}^{n+1} \right)_{ij}\cdot \vect{S}_{ij}\\
\\
\text{\ \ \ \ } - \displaystyle\sum\limits_{j\in V(i)}
\left(K^n \gradv\left(\frac{1}{2}(u^2)^{n+1}
+ \varepsilon_{sup}^{n+\frac{1}{2}}\right)\right)_{ij}\cdot\vect{S}_{ij}
+ \Omega_i\rho_i^{n+1}{\Phi_v}_i\\
\end{array}
\end{equation}


%.................................
\subsubsection*{Discrétisation de la partie ``convective''}
%.................................

La valeur à la face s'écrit~:
\begin{equation}
\left(e^{n+1} \vect{Q}_{ac}^{n+1}\right)_{ij} \cdot \vect{S}_{ij}
= e_{ij}^{n+1}(\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij}
\end{equation}
avec un décentrement sur la valeur de $e^{n+1}$ aux faces~:
\begin{equation}
\begin{array}{lllll}
e_{ij}^{n+1}
& = & e_i^{n+1}
& \text{si\ } & (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} \geqslant 0 \\
& = & e_j^{n+1}
& \text{si\ } & (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} < 0 \\
\end{array}
\end{equation}
que l'on peut noter~:
\begin{equation}
 e_{ij}^{n+1}
 = \beta_{ij}e_i^{n+1} + (1-\beta_{ij})e_j^{n+1}
\end{equation}
avec
\begin{equation}
\left\{\begin{array}{lll}
\beta_{ij} = 1 & \text{si\ }
& (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} \geqslant 0 \\
\beta_{ij} = 0 & \text{si\ }
& (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} < 0 \\
\end{array}\right.
\end{equation}


%.................................
\subsubsection*{Discrétisation de la partie ``diffusive''}
%.................................

La valeur à la face s'écrit~:
\begin{equation}
\begin{array}{c}
\left(K^n\gradv(e^{n+1})\right)_{ij}\cdot\vect{S}_{ij}
= K_{ij}^n
\displaystyle \left( \frac{\partial e}{\partial n} \right)^{n+1}_{ij}S_{ij}\\
\text{et}\\
\left(K^n \gradv\left(\frac{1}{2}(u^2)^{n+1}
+ \varepsilon_{sup}^{n+\frac{1}{2}}\right)\right)_{ij}\cdot\vect{S}_{ij}
= K_{ij}^n
\displaystyle \left( \frac{\partial \left(\frac{1}{2} u^2
+ \varepsilon_{sup}\right)}{\partial n} \right)^{n+\frac{1}{2}}_{ij}S_{ij}
\end{array}
\end{equation}
avec une interpolation linéaire pour
$K^n$ aux faces (et en pratique, $\alpha_{ij}=\frac{1}{2}$)~:
\begin{equation}
K_{ij}^n
= \alpha_{ij}K_{i}^n+(1-\alpha_{ij})K_{j}^n
\end{equation}
et un schéma centré avec reconstruction pour le gradient normal aux faces~:
\begin{equation}
\displaystyle \left( \frac{\partial e}{\partial n} \right)^{n+1}_{ij}
= \displaystyle\frac{e_{J'}^{n+1} - e_{I'}^{n+1}}{\overline{I'J'}}
\quad \text{et} \quad
\displaystyle \left( \frac{\partial \left(\frac{1}{2} u^2
+ \varepsilon_{sup}\right)}{\partial n} \right)^{n+\frac{1}{2}}_{ij}
= \displaystyle\frac{(\frac{1}{2} u^2
+ \varepsilon_{sup})_{J'}^{n+\frac{1}{2}} - (\frac{1}{2} u^2
+ \varepsilon_{sup})_{I'}^{n+\frac{1}{2}}}{\overline{I'J'}}
\end{equation}



%.................................
\subsubsection*{Discrétisation de la puissance des forces de pression}
%.................................

Ce terme
est issu du terme convectif, on le discrétise donc de la même façon.

\begin{equation}
\left(\displaystyle\frac{\widetilde{P}}{\rho^{n+1}}\
\vect{Q}_{ac}^{n+1}\right)_{ij} \cdot \vect{S}_{ij}
= \left(\displaystyle\frac{\widetilde{P}}{\rho^{n+1}}\right)_{ij}
(\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij}
\end{equation}

avec un décentrement sur la valeur de
$\displaystyle\frac{P}{\rho}$ aux faces~:
\begin{equation}
\begin{array}{lll}
\left(\displaystyle\frac{\widetilde{P}}{\rho^{n+1}}\right)_{ij}
 = \beta_{ij}\displaystyle\frac{\widetilde{P}_i}{\rho^{n+1}_i}
+ (1-\beta_{ij})\displaystyle\frac{\widetilde{P}_j}{\rho^{n+1}_j}
& \text{avec}
& \left\{\begin{array}{lll}
\beta_{ij} = 1 & \text{si\ }
& (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} \geqslant 0 \\
\beta_{ij} = 0 & \text{si\ }
& (\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij} < 0 \\
\end{array}\right.
\end{array}
\end{equation}



%.................................
\subsubsection*{Discrétisation de la puissance des forces visqueuses}
%.................................

On calcule les termes dans les cellules puis on utilise une
interpolation linéaire (on utilise
$\alpha_{ij}=\frac{1}{2}$ dans la relation ci-dessous)~:
\begin{equation}
\left((\tens{\Sigma}^v)^{n+1} \vect{u}^{n+1} \right)_{ij}\cdot \vect{S}_{ij}
= \left\{\alpha_{ij} \left((\tens{\Sigma}^v)^{n+1} \vect{u}^{n+1}\right)_i
+ (1-\alpha_{ij}) \left((\tens{\Sigma}^v)^{n+1} \vect{u}^{n+1}\right)_j
\right\} \cdot \vect{S}_{ij}
\end{equation}


%.................................
\subsubsection*{Remarques}
%.................................


Les termes ``convectifs'' associés à
$\displaystyle\dive\left(\left(e^{n+1}+\frac{\widetilde{P}}{\rho^{n+1}}\right)\,
\vect{Q}_{ac}^{n+1}\right)$ sont calculés avec un décentrement amont
(consistant, d'ordre 1 en espace). Les valeurs utilisées sont bien prises au
centre de la cellule amont ($e_i$, $P_i$, $\rho_i$) et non pas au projeté $I'$
du centre de la cellule sur la normale à la face passant par son centre de
gravité (sur un cas test en triangles, l'utilisation de $P_I'$ et de $\rho_I'$
pour le terme de transport de pression a conduit à un résultat
insatisfaisant, mais des corrections ont été apportées aux sources depuis
et il serait utile de vérifier que cette conclusion n'est pas remise en question).

Les termes diffusifs associés à
$\displaystyle\dive\left(K\,\grad\left(e+\frac{1}{2} u^2 +
\varepsilon_{sup}\right)\right)$ sont calculés en utilisant des valeurs aux
faces reconstruites pour s'assurer de la consistance du schéma.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Après une étape de gestion de la mémoire (\fort{memcfe}), on calcule les
différents termes sources (au centre des cellules)~:
\begin{itemize}
\item source volumique de chaleur (\fort{cs\_user\_source\_terms}),
\item source associée aux sources de masse (\fort{catsma}),
\item source associée à l'accumulation de masse $\dive\,\vect{Q}_{ac}$ (directement dans \fort{cfener}),
\item dissipation visqueuse (\fort{cfdivs}),
\item transport de pression (directement dans \fort{cfener}),
\item puissance de la pesanteur (directement dans \fort{cfener}),
\item termes diffusifs en $\displaystyle\dive\left(K\,\grad\left(\frac{1}{2} u^2 +
\varepsilon_{sup}\right)\right)$ (calcul de $\varepsilon_{sup}$ par
\fort{uscfth}, puis calcul du terme diffusif directement dans \fort{cfener}).
\end{itemize}

\bigskip
Le système (\ref{Cfbl_Cfener_eq_energie_totale_discrete_cfener}) est résolu par une méthode
d'incrément et résidu  en utilisant une méthode de Jacobi (\fort{cfcdts}).

L'impression des bornes et
la limitation éventuelle de l'énergie sont ensuite effectuées par
\fort{clpsca} suivi de \fort{uscfth} (intervention utilisateur optionnelle).

On actualise enfin la pression et on calcule la
température (\fort{uscfth}).

Pour
terminer, en parallèle ou en périodique, on échange les variables
pression, énergie et température.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Points à traiter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% propose en patch 1.2.1

%Corriger \fort{cfener} dans lequel \var{W1} produit par \fort{uscfth} est
%écrasé par \fort{grdcel}, causant probablement des dégâts
%dans les cas où le gradient de l'énergie cinétique dans la direction $x$
%est sensiblement non nul sur des faces de bord dont la normale a une
%composante en $x$ et lorsque la conductivité n'est pas négligeable.


\etape{Choix de $\widetilde{P}$}
En standard, on utilise $\widetilde{P}=P^n$, mais ce n'est pas le seul choix
possible. On pourrait étudier le comportement de l'algorithme avec $P^{Pred}$ et
$P^{n+1}$ (avec $P^{n+1}$, en particulier,
$\displaystyle\frac{\widetilde{P}}{\rho^{n+1}}$
est évalué avec la masse volumique et l'énergie prises au même instant).

\etape{Terme source dans l'équation de l'énergie}
La présence d'un terme source externe dans l'équation de l'énergie génère des
oscillations de vitesse qu'il est important d'analyser et de comprendre.
